Podsumowanie analizy

Celem analizy było określenie głównych przyczyn stopniowego zmniejszania się długości śledzi oceanicznych wyławianych w Europie. Analiza doprowadziła do wniosków, że istotny wpływ na długość miała temperatura przy powierzchni wody i to ona mogła doprowadzić do zmian długości. Przy wzroscie temperatury z roku na rok długość śledzia stopniowo malała. Duży wpływ na długość śledzia miały również wartości związane z dostępnością planktonu, co jest logiczne gdyż plankton jest pożywieniem dla śledzi. Wszystkie te wartości mogą jednak być związane z oscylacją północnoatlantycką, która jest zjawiskiem meteorologicznym i może mieć związek z temperaturą wody przy powierzchni i poprzez cyrkulację wody oceanicznej z dostępnością planktonu. Wartość ta jest przeciwnie skorelowana z długością, jednakże żaden regresor nie uznał jej za istotną. W danych brakuje informacji o roku pomiaru co utrudnia analizę, jednakże można ją wyłuskać z atrybutów rocznego narybku i rocznego natężenia połowów w rejonie.

Wprowadzenie

Opis problemu

Celem projektu jest określenie głównych przyczyn stopniowego zmniejszania się długości śledzi oceanicznych wyławionych w Europie.

Opis danych

Zbiór danych do analizy składa się z pomiarów śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Opis zmiennych:

Wszystkie dane, są danymi liczbowymi, nie występują żadne dane kategoryczne.

length: długość złowionego śledzia [cm];
cfin1: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1];
cfin2: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2];
chel1: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1];
chel2: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2];
lcop1: dostępność planktonu [zagęszczenie widłonogów gat. 1];
lcop2: dostępność planktonu [zagęszczenie widłonogów gat. 2];
fbar: natężenie połowów w regionie [ułamek pozostawionego narybku];
recr: roczny narybek [liczba śledzi];
cumf: łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku];
totaln: łączna liczba ryb złowionych w ramach połowu [liczba śledzi];
sst: temperatura przy powierzchni wody [°C];
sal: poziom zasolenia wody [Knudsen ppt];
xmonth: miesiąc połowu [numer miesiąca];
nao: oscylacja północnoatlantycka [mb].

Inicjalizacja środowiska

library(dplyr)
library(ggplot2)
library(GGally) # ggpairplot
library(DT)
library(mice) # for missing data
library(ggpubr) # ggarange
library(corrplot) # ggarange
library(plotly)
library(caret)
library(gbm)
library(doMC)

columns_names <- c("length" = "długość złowionego śledzia [cm]",
  "cfin1" = "dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]",
  "cfin2" = "dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]",
  "chel1" = "dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]",
  "chel2" = "dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]",
  "lcop1" = "dostępność planktonu [zagęszczenie widłonogów gat. 1]",
  "lcop2" = "dostępność planktonu [zagęszczenie widłonogów gat. 2]",
  "fbar" = "natężenie połowów w regionie [ułamek pozostawionego narybku]",
  "recr" = "roczny narybek [liczba śledzi]",
  "cumf" = "łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]",
  "totaln" = "łączna liczba ryb złowionych w ramach połowu [liczba śledzi]",
  "sst" = "temperatura przy powierzchni wody [°C]",
  "sal"  = "poziom zasolenia wody [Knudsen ppt]",
  "xmonth" = "miesiąc połowu [numer miesiąca]",
  "nao" = "oscylacja północnoatlantycka [mb]" )

prettyTable <- function(table_df, round_columns=numeric(), round_digits=2) {
    DT::datatable(table_df, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print'))) %>%
    formatRound(round_columns, round_digits)
}

Analiza

Wczytanie i wstepne oczyszczenie danych

Pierwszym krokiem jest wczytanie danych oraz bliższe przyjżenie się im. Ponieważ w pliku csv brakujące rekordy zapisane są jako “?” trzeba je poprawnie wczytać. Brakujące dane występują w kolumnach cfin1, cfin2, chel1, chel2, lcop1, lcop2 i lsst. Jako że kolumna x to po prostu numer pomiaru, nie będzie nam potrzebna.

set.seed(23)
data <- read.csv("sledzie.csv", na.strings = c("?"))
lapply(data, function(x) any(is.na(x)))

## $X
## [1] FALSE
## 
## $length
## [1] FALSE
## 
## $cfin1
## [1] TRUE
## 
## $cfin2
## [1] TRUE
## 
## $chel1
## [1] TRUE
## 
## $chel2
## [1] TRUE
## 
## $lcop1
## [1] TRUE
## 
## $lcop2
## [1] TRUE
## 
## $fbar
## [1] FALSE
## 
## $recr
## [1] FALSE
## 
## $cumf
## [1] FALSE
## 
## $totaln
## [1] FALSE
## 
## $sst
## [1] TRUE
## 
## $sal
## [1] FALSE
## 
## $xmonth
## [1] FALSE
## 
## $nao
## [1] FALSE

data <- subset(data, select = -c(X))

summary(data)

##      length         cfin1             cfin2             chel1       
##  Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778   1st Qu.: 2.469  
##  Median :25.5   Median : 0.1111   Median : 0.7012   Median : 5.750  
##  Mean   :25.3   Mean   : 0.4458   Mean   : 2.0248   Mean   :10.006  
##  3rd Qu.:26.5   3rd Qu.: 0.3333   3rd Qu.: 1.7936   3rd Qu.:11.500  
##  Max.   :32.5   Max.   :37.6667   Max.   :19.3958   Max.   :75.000  
##                 NA's   :1581      NA's   :1536      NA's   :1555    
##      chel2            lcop1              lcop2             fbar       
##  Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849   Min.   :0.0680  
##  1st Qu.:13.427   1st Qu.:  2.5479   1st Qu.:17.808   1st Qu.:0.2270  
##  Median :21.673   Median :  7.0000   Median :24.859   Median :0.3320  
##  Mean   :21.221   Mean   : 12.8108   Mean   :28.419   Mean   :0.3304  
##  3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232   3rd Qu.:0.4560  
##  Max.   :57.706   Max.   :115.5833   Max.   :68.736   Max.   :0.8490  
##  NA's   :1556     NA's   :1653       NA's   :1591                     
##       recr              cumf             totaln             sst       
##  Min.   : 140515   Min.   :0.06833   Min.   : 144137   Min.   :12.77  
##  1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068   1st Qu.:13.60  
##  Median : 421391   Median :0.23191   Median : 539558   Median :13.86  
##  Mean   : 520367   Mean   :0.22981   Mean   : 514973   Mean   :13.87  
##  3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351   3rd Qu.:14.16  
##  Max.   :1565890   Max.   :0.39801   Max.   :1015595   Max.   :14.73  
##                                                        NA's   :1584   
##       sal            xmonth            nao          
##  Min.   :35.40   Min.   : 1.000   Min.   :-4.89000  
##  1st Qu.:35.51   1st Qu.: 5.000   1st Qu.:-1.89000  
##  Median :35.51   Median : 8.000   Median : 0.20000  
##  Mean   :35.51   Mean   : 7.258   Mean   :-0.09236  
##  3rd Qu.:35.52   3rd Qu.: 9.000   3rd Qu.: 1.63000  
##  Max.   :35.61   Max.   :12.000   Max.   : 5.08000  
##

paste0("Ilość badanych przypadków: ", nrow(data), ", ilość zmiennych: ", ncol(data))

## [1] "Ilosc badanych przypadków: 52582, ilosc zmiennych: 15"

head(data, 10)

##    length   cfin1   cfin2   chel1    chel2   lcop1    lcop2  fbar   recr
## 1    23.0 0.02778 0.27785 2.46875       NA 2.54787 26.35881 0.356 482831
## 2    22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3    25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4    25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5    24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 6    22.0 0.02778 0.27785 2.46875 21.43548 2.54787       NA 0.356 482831
## 7    24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 8    23.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 9    22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 10   22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
##         cumf   totaln      sst      sal xmonth nao
## 1  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 2  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 3  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 4  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 5  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 6  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 7  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 8  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 9  0.3059879 267380.8 14.30693 35.51234      7 2.8
## 10 0.3059879 267380.8 14.30693 35.51234      7 2.8

W celu poradzenia sobie z brakującymi danymi zostaje wykorzystana metoda MICE (Multivariate Imputation via Chained Equations). Metoda ta pozwala na pozbycie się brakujących danych z utrzymaniem podobnego rozkładu.

plot_histograms <- function(data, columns=NULL, ncol=1, nrow=1, createName=function(x) x, bins=30){
  if(is.null(columns)){
    columns <- names(data)
  } 
  
  i <- 1
  plots <- list()
  for (column_name in columns){
    p <- ggplot(data, aes_string(x = column_name)) + 
      geom_histogram(color="black", fill="orange", bins=bins) +
      # geom_density() +
      ggtitle(createName(column_name)) + 
      ylab("Liczba obserwacji")
    plots[[i]] <- p
    i <- i + 1
  }
  figure <- ggarrange(plotlist = plots, ncol = ncol, nrow = nrow)
  figure
}


  

missing_data_columns <- c("cfin1", "cfin2", "chel1", "chel2", "lcop1", "lcop2", "sst")
# plot_histograms(data, columns = missing_data_columns, createName = function(x) paste0("Histogram ", x, " przed pozbyciem się brakujących danych"), nrow=2)
# mice_plot <- aggr(data, col=c('navyblue','yellow'),
#                     numbers=TRUE, sortVars=TRUE,
#                     labels=names(data), cex.axis=.7,
#                     gap=3, ylab=c("Missing data","Pattern"))
completed_data <- data %>%
  mice(m=3,  method = 'cart', seed = 63) %>%
  complete(2)

## 
##  iter imp variable
##   1   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   1   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   1   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst

# plot_histograms(completed_data, columns = missing_data_columns, createName = function(x) paste0("Histogram ", x, " po pozbyciu się brakujących danych"), nrow=2)

Eksploracja danych

Podsumowanie danych ich podstawowych statystyk.

summary(completed_data)

##      length         cfin1             cfin2             chel1       
##  Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778   1st Qu.: 2.469  
##  Median :25.5   Median : 0.1111   Median : 0.7012   Median : 5.750  
##  Mean   :25.3   Mean   : 0.4460   Mean   : 2.0255   Mean   :10.000  
##  3rd Qu.:26.5   3rd Qu.: 0.3333   3rd Qu.: 1.7936   3rd Qu.:11.500  
##  Max.   :32.5   Max.   :37.6667   Max.   :19.3958   Max.   :75.000  
##      chel2            lcop1              lcop2             fbar       
##  Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849   Min.   :0.0680  
##  1st Qu.:13.427   1st Qu.:  2.5479   1st Qu.:17.808   1st Qu.:0.2270  
##  Median :21.673   Median :  7.0000   Median :24.859   Median :0.3320  
##  Mean   :21.219   Mean   : 12.8060   Mean   :28.422   Mean   :0.3304  
##  3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232   3rd Qu.:0.4560  
##  Max.   :57.706   Max.   :115.5833   Max.   :68.736   Max.   :0.8490  
##       recr              cumf             totaln             sst       
##  Min.   : 140515   Min.   :0.06833   Min.   : 144137   Min.   :12.77  
##  1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068   1st Qu.:13.60  
##  Median : 421391   Median :0.23191   Median : 539558   Median :13.86  
##  Mean   : 520367   Mean   :0.22981   Mean   : 514973   Mean   :13.87  
##  3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351   3rd Qu.:14.16  
##  Max.   :1565890   Max.   :0.39801   Max.   :1015595   Max.   :14.73  
##       sal            xmonth            nao          
##  Min.   :35.40   Min.   : 1.000   Min.   :-4.89000  
##  1st Qu.:35.51   1st Qu.: 5.000   1st Qu.:-1.89000  
##  Median :35.51   Median : 8.000   Median : 0.20000  
##  Mean   :35.51   Mean   : 7.258   Mean   :-0.09236  
##  3rd Qu.:35.52   3rd Qu.: 9.000   3rd Qu.: 1.63000  
##  Max.   :35.61   Max.   :12.000   Max.   : 5.08000

Histogramy wartości atrybutów w zbiorze danych.

Jak można zauważyć wartości długości śledzia przyjmują rozkład normalny ze średnią na poziomie 25.3 cm i przyjmują wartości w granicach 19-32.5cm.

Najwięcej połowów występuje w okresie letnim (czerwiec - wrzesień) oraz w październiku.

Pozostałe atrybuty zostały opisane na poniższych histogramach.

plot_histograms(completed_data, createName = function(x) paste0("Histogram - ", columns_names[x]), nrow=1, bins=40)

## $`1`

## 
## $`2`

## 
## $`3`

## 
## $`4`

## 
## $`5`

## 
## $`6`

## 
## $`7`

## 
## $`8`

## 
## $`9`

## 
## $`10`

## 
## $`11`

## 
## $`12`

## 
## $`13`

## 
## $`14`

## 
## $`15`

## 
## attr(,"class")
## [1] "list"      "ggarrange"

Poniżej została przedstawiona korelacja pomiędzy atrybutami. Można zauważyć, że występuje wysoka korelacja pomiędzy parami atrybutów chel1-lcop1 oraz chel2-lcop2, ponieważ związane są one z dostępnością planktonu. Korelacja występuje także pomiędzy fbar-cumf, ponieważ obydwa atrybuty związane są z natężeniem popłowu w regionie. Można te pary zredukować, aby przeciwdziałać klątwie wielowymiarowości - jednakże tutaj zostaną one zostawione bez zmian.

corrplot(cor(completed_data), type = "upper",  tl.col = "black", tl.srt = 45, method="number")

ggplot(completed_data, aes(x=chel1, y=lcop1)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  ggtitle("Zależność pomiędzy atrybutem chel1 i lcop1")

linear_model_1 <- lm(chel1 ~ lcop1, data = completed_data)


ggplot(completed_data, aes(x=chel2, y=lcop2)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  ggtitle("Zależność pomiędzy atrybutem chel2 i lcop2")

linear_model_2 <- lm(chel2 ~ lcop2, data = completed_data)


ggplot(completed_data, aes(x=fbar, y=cumf)) + 
  geom_point() + 
  geom_smooth(method=lm) +
  ggtitle("Zależność pomiędzy atrybutem fbar i cumf")

linear_model_3 <- lm(fbar ~ cumf, data = completed_data)

Poniżej przedstawiona została analiza długości śledzia na przestrzeni czasu. Jednakże czas nie był przedstawiony jawnie, ale wiadomo, że dane są ustawione chronologicznie dzięki czemu, zbiór można podzielić i na każdym z podziałów policzyć średnią. Można zauważyć, że na początku długość śledzia na przestrzeni czasu wzrastała od średniej wielkości w oklicach 24.5 cm do wielkości wielkości przekraczającej 27cm. Jednakże później trend się odwrócił i ostatecznie na sam koniec średnia ta była poniżej 24 cm.

aggregate_data <- function(data, cases_per_aggregate){
  new_aggregated_data <- aggregate(data, 
                        list(rep(1:(nrow(data)%/%cases_per_aggregate+1), each=cases_per_aggregate, len=nrow(data))),
                        mean) # Try using mean without outliners
  new_aggregated_data$time_period = c(1:nrow(new_aggregated_data))
  new_aggregated_data
}


aggregated_data <- aggregate_data(completed_data, 200)
ggplot(aggregated_data, aes(time_period, length)) +
  geom_point() + 
  geom_smooth() +
  ggtitle("Wykres zmiany długości śledzia na przestrzeni czasu") +
  ylab("Długość śledzia") +
  xlab("Przedział czasowy")

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot_ly(data = aggregated_data, x = ~time_period, y = ~length,
            text = ~paste("Length: ", length))

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

# %>% layout(
#               scene = list(
#                  xaxis = list(title = "Przedział czasowy"),
#                  yaxis = list(title = "Długość śledzia")
#               )
#             )

Regresor

train_index <- createDataPartition(completed_data$length, p = .7, 
                                  list = FALSE, 
                                  times = 1)
training <- completed_data[ train_index,]
testing  <- completed_data[-train_index,]

r_squared <- function (x, y) cor(x, y) ^ 2

registerDoMC(cores=8)

ctrl <- trainControl(
  method = "repeatedcv",
  number = 5,
  repeats = 3
  )

evaluate_model <- function(model, test_data) {
  print(model)
  plot <- ggplot(varImp(model))
  print(varImp(model))
  print(plot)
  predicted <- predict(model, newdata = test_data)
  rmse_ <- RMSE(predicted, test_data$length)
  rsquared <- r_squared(predicted, test_data$length)
  print(paste0("RMSE na zbiorze testowym: ", rmse_, " R-Squared na zbiorze testowym ", rsquared))
  c(rmse_, rsquared)
}

Elastic Net

modelLookup(model='glmnet')

##    model parameter                    label forReg forClass probModel
## 1 glmnet     alpha        Mixing Percentage   TRUE     TRUE      TRUE
## 2 glmnet    lambda Regularization Parameter   TRUE     TRUE      TRUE

glmnet_grid = expand.grid(
  alpha = 0:1,
  lambda = seq(0.0001, 1, length = 100)
)

set.seed(23)

glmnet <- train(length ~ .,
             data = training,
             method = "glmnet",
             preProc = c("center", "scale"),
             metric = "RMSE",
             trControl = ctrl,
             tuneGrid = glmnet_grid,
             importance=TRUE)

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info =
## trainInfo, : There were missing values in resampled performance measures.

glmnet_metrics = evaluate_model(glmnet, testing)

## glmnet 
## 
## 36810 samples
##    14 predictor
## 
## Pre-processing: centered (14), scaled (14) 
## Resampling: Cross-Validated (5 fold, repeated 3 times) 
## Summary of sample sizes: 29447, 29449, 29449, 29447, 29448, 29448, ... 
## Resampling results across tuning parameters:
## 
##   alpha  lambda  RMSE      Rsquared   MAE     
##   0      0.0001  1.370182  0.3143406  1.091332
##   0      0.0102  1.370182  0.3143406  1.091332
##   0      0.0203  1.370182  0.3143406  1.091332
##   0      0.0304  1.370182  0.3143406  1.091332
##   0      0.0405  1.370182  0.3143406  1.091332
##   0      0.0506  1.370182  0.3143406  1.091332
##   0      0.0607  1.370182  0.3143406  1.091332
##   0      0.0708  1.370182  0.3143406  1.091332
##   0      0.0809  1.370975  0.3137485  1.091997
##   0      0.0910  1.372204  0.3128238  1.093068
##   0      0.1011  1.373432  0.3118982  1.094216
##   0      0.1112  1.374655  0.3109730  1.095383
##   0      0.1213  1.375870  0.3100524  1.096521
##   0      0.1314  1.377076  0.3091374  1.097612
##   0      0.1415  1.378265  0.3082344  1.098650
##   0      0.1516  1.379429  0.3073522  1.099638
##   0      0.1617  1.380579  0.3064784  1.100609
##   0      0.1718  1.381721  0.3056095  1.101582
##   0      0.1819  1.382818  0.3047808  1.102528
##   0      0.1920  1.383915  0.3039490  1.103488
##   0      0.2021  1.384983  0.3031414  1.104450
##   0      0.2122  1.386037  0.3023462  1.105412
##   0      0.2223  1.387075  0.3015634  1.106367
##   0      0.2324  1.388090  0.3008004  1.107301
##   0      0.2425  1.389083  0.3000564  1.108214
##   0      0.2526  1.390069  0.2993171  1.109116
##   0      0.2627  1.391016  0.2986134  1.109975
##   0      0.2728  1.391980  0.2978926  1.110840
##   0      0.2829  1.392891  0.2972216  1.111651
##   0      0.2930  1.393809  0.2965424  1.112462
##   0      0.3031  1.394718  0.2958696  1.113262
##   0      0.3132  1.395580  0.2952419  1.114014
##   0      0.3233  1.396456  0.2945998  1.114773
##   0      0.3334  1.397314  0.2939733  1.115509
##   0      0.3435  1.398135  0.2933822  1.116206
##   0      0.3536  1.398968  0.2927788  1.116905
##   0      0.3637  1.399796  0.2921783  1.117598
##   0      0.3738  1.400574  0.2916267  1.118248
##   0      0.3839  1.401362  0.2910647  1.118907
##   0      0.3940  1.402160  0.2904927  1.119576
##   0      0.4041  1.402916  0.2899607  1.120212
##   0      0.4142  1.403659  0.2894412  1.120834
##   0      0.4243  1.404410  0.2889128  1.121458
##   0      0.4344  1.405166  0.2883789  1.122081
##   0      0.4445  1.405873  0.2878921  1.122659
##   0      0.4546  1.406576  0.2874098  1.123228
##   0      0.4647  1.407286  0.2869200  1.123797
##   0      0.4748  1.408002  0.2864233  1.124366
##   0      0.4849  1.408681  0.2859627  1.124903
##   0      0.4950  1.409343  0.2855180  1.125424
##   0      0.5051  1.410012  0.2850668  1.125947
##   0      0.5152  1.410687  0.2846092  1.126474
##   0      0.5253  1.411360  0.2841532  1.127002
##   0      0.5354  1.411988  0.2837410  1.127506
##   0      0.5455  1.412617  0.2833277  1.128024
##   0      0.5556  1.413252  0.2829089  1.128557
##   0      0.5657  1.413892  0.2824845  1.129103
##   0      0.5758  1.414529  0.2820619  1.129652
##   0      0.5859  1.415124  0.2816807  1.130171
##   0      0.5960  1.415716  0.2813025  1.130691
##   0      0.6061  1.416312  0.2809196  1.131218
##   0      0.6162  1.416913  0.2805321  1.131752
##   0      0.6263  1.417518  0.2801402  1.132292
##   0      0.6364  1.418103  0.2797678  1.132816
##   0      0.6465  1.418660  0.2794233  1.133318
##   0      0.6566  1.419220  0.2790749  1.133823
##   0      0.6667  1.419784  0.2787226  1.134330
##   0      0.6768  1.420352  0.2783663  1.134839
##   0      0.6869  1.420923  0.2780063  1.135349
##   0      0.6970  1.421481  0.2776595  1.135845
##   0      0.7071  1.422006  0.2773437  1.136314
##   0      0.7172  1.422535  0.2770248  1.136783
##   0      0.7273  1.423068  0.2767024  1.137252
##   0      0.7374  1.423603  0.2763767  1.137722
##   0      0.7475  1.424142  0.2760476  1.138192
##   0      0.7576  1.424682  0.2757175  1.138660
##   0      0.7677  1.425196  0.2754119  1.139106
##   0      0.7778  1.425691  0.2751242  1.139535
##   0      0.7879  1.426190  0.2748338  1.139964
##   0      0.7980  1.426691  0.2745405  1.140394
##   0      0.8081  1.427195  0.2742443  1.140825
##   0      0.8182  1.427702  0.2739453  1.141257
##   0      0.8283  1.428211  0.2736438  1.141689
##   0      0.8384  1.428708  0.2733533  1.142112
##   0      0.8485  1.429176  0.2730912  1.142511
##   0      0.8586  1.429646  0.2728278  1.142911
##   0      0.8687  1.430117  0.2725620  1.143313
##   0      0.8788  1.430591  0.2722939  1.143717
##   0      0.8889  1.431068  0.2720233  1.144122
##   0      0.8990  1.431546  0.2717502  1.144528
##   0      0.9091  1.432027  0.2714752  1.144935
##   0      0.9192  1.432498  0.2712089  1.145333
##   0      0.9293  1.432941  0.2709697  1.145707
##   0      0.9394  1.433383  0.2707307  1.146079
##   0      0.9495  1.433828  0.2704896  1.146452
##   0      0.9596  1.434274  0.2702466  1.146825
##   0      0.9697  1.434723  0.2700015  1.147198
##   0      0.9798  1.435173  0.2697543  1.147572
##   0      0.9899  1.435626  0.2695050  1.147946
##   0      1.0000  1.436079  0.2692550  1.148319
##   1      0.0001  1.363552  0.3192588  1.084409
##   1      0.0102  1.365472  0.3177033  1.085918
##   1      0.0203  1.369221  0.3148282  1.089374
##   1      0.0304  1.374165  0.3110152  1.093505
##   1      0.0405  1.379621  0.3068474  1.099438
##   1      0.0506  1.386537  0.3011004  1.106289
##   1      0.0607  1.392899  0.2959747  1.112238
##   1      0.0708  1.399495  0.2905558  1.118008
##   1      0.0809  1.407027  0.2838729  1.124723
##   1      0.0910  1.415296  0.2760606  1.132537
##   1      0.1011  1.423924  0.2675317  1.140801
##   1      0.1112  1.433152  0.2578987  1.149419
##   1      0.1213  1.441206  0.2494461  1.156718
##   1      0.1314  1.446138  0.2449423  1.161262
##   1      0.1415  1.448712  0.2436123  1.163650
##   1      0.1516  1.451009  0.2427450  1.165754
##   1      0.1617  1.453433  0.2417993  1.167946
##   1      0.1718  1.456002  0.2407395  1.170163
##   1      0.1819  1.458721  0.2395443  1.172424
##   1      0.1920  1.461590  0.2382005  1.174756
##   1      0.2021  1.464608  0.2366931  1.177180
##   1      0.2122  1.467774  0.2350070  1.180064
##   1      0.2223  1.471086  0.2331238  1.183065
##   1      0.2324  1.474546  0.2310234  1.186110
##   1      0.2425  1.478145  0.2286967  1.189240
##   1      0.2526  1.481634  0.2265382  1.192206
##   1      0.2627  1.485212  0.2242264  1.195414
##   1      0.2728  1.488566  0.2223093  1.198411
##   1      0.2829  1.491748  0.2207225  1.201208
##   1      0.2930  1.495038  0.2189623  1.204101
##   1      0.3031  1.498436  0.2170113  1.207070
##   1      0.3132  1.501941  0.2148499  1.210146
##   1      0.3233  1.505550  0.2124620  1.213486
##   1      0.3334  1.509226  0.2098939  1.216839
##   1      0.3435  1.513005  0.2070578  1.220220
##   1      0.3536  1.516780  0.2041574  1.223578
##   1      0.3637  1.519431  0.2036212  1.225840
##   1      0.3738  1.522063  0.2032341  1.228081
##   1      0.3839  1.524757  0.2028211  1.230388
##   1      0.3940  1.527385  0.2026949  1.232609
##   1      0.4041  1.530021  0.2026949  1.234872
##   1      0.4142  1.532719  0.2026949  1.237169
##   1      0.4243  1.535479  0.2026949  1.239466
##   1      0.4344  1.538300  0.2026949  1.241771
##   1      0.4445  1.541182  0.2026949  1.244114
##   1      0.4546  1.544125  0.2026949  1.246473
##   1      0.4647  1.547128  0.2026949  1.248842
##   1      0.4748  1.550191  0.2026949  1.251314
##   1      0.4849  1.553314  0.2026949  1.253911
##   1      0.4950  1.556496  0.2026949  1.256537
##   1      0.5051  1.559737  0.2026949  1.259215
##   1      0.5152  1.563036  0.2026949  1.262035
##   1      0.5253  1.566394  0.2026949  1.264921
##   1      0.5354  1.569809  0.2026949  1.267861
##   1      0.5455  1.573282  0.2026949  1.270836
##   1      0.5556  1.576812  0.2026949  1.273878
##   1      0.5657  1.580399  0.2026949  1.276973
##   1      0.5758  1.584042  0.2026949  1.280119
##   1      0.5859  1.587741  0.2026949  1.283340
##   1      0.5960  1.591495  0.2026949  1.286631
##   1      0.6061  1.595305  0.2026949  1.290000
##   1      0.6162  1.599169  0.2026949  1.293437
##   1      0.6263  1.603087  0.2026949  1.296960
##   1      0.6364  1.607060  0.2026949  1.300528
##   1      0.6465  1.611086  0.2026949  1.304142
##   1      0.6566  1.615165  0.2026949  1.307781
##   1      0.6667  1.619297  0.2026949  1.311441
##   1      0.6768  1.623481  0.2026949  1.315195
##   1      0.6869  1.627717  0.2026949  1.318972
##   1      0.6970  1.632005  0.2026949  1.322749
##   1      0.7071  1.636343  0.2026949  1.326527
##   1      0.7172  1.640733  0.2026949  1.330304
##   1      0.7273  1.645172  0.2026949  1.334081
##   1      0.7374  1.649614  0.2015148  1.337818
##   1      0.7475  1.652483  0.1944352  1.340208
##   1      0.7576  1.652485        NaN  1.340209
##   1      0.7677  1.652485        NaN  1.340209
##   1      0.7778  1.652485        NaN  1.340209
##   1      0.7879  1.652485        NaN  1.340209
##   1      0.7980  1.652485        NaN  1.340209
##   1      0.8081  1.652485        NaN  1.340209
##   1      0.8182  1.652485        NaN  1.340209
##   1      0.8283  1.652485        NaN  1.340209
##   1      0.8384  1.652485        NaN  1.340209
##   1      0.8485  1.652485        NaN  1.340209
##   1      0.8586  1.652485        NaN  1.340209
##   1      0.8687  1.652485        NaN  1.340209
##   1      0.8788  1.652485        NaN  1.340209
##   1      0.8889  1.652485        NaN  1.340209
##   1      0.8990  1.652485        NaN  1.340209
##   1      0.9091  1.652485        NaN  1.340209
##   1      0.9192  1.652485        NaN  1.340209
##   1      0.9293  1.652485        NaN  1.340209
##   1      0.9394  1.652485        NaN  1.340209
##   1      0.9495  1.652485        NaN  1.340209
##   1      0.9596  1.652485        NaN  1.340209
##   1      0.9697  1.652485        NaN  1.340209
##   1      0.9798  1.652485        NaN  1.340209
##   1      0.9899  1.652485        NaN  1.340209
##   1      1.0000  1.652485        NaN  1.340209
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 1e-04.
## glmnet variable importance
## 
##        Overall
## fbar   100.000
## cumf    95.433
## sst     54.392
## lcop1   30.601
## chel1   12.388
## totaln  11.977
## lcop2   11.788
## cfin1   11.436
## recr     9.039
## cfin2    6.426
## nao      4.945
## chel2    1.719
## xmonth   1.099
## sal      0.000

## [1] "RMSE na zbiorze testowym: 1.35944937812852 R-Squared na zbiorze testowym 0.324293505084823"

Random Forest

set.seed(23)
modelLookup(model='rf')

##   model parameter                         label forReg forClass probModel
## 1    rf      mtry #Randomly Selected Predictors   TRUE     TRUE      TRUE

rf_grid <- expand.grid(mtry=c(2, 3, 4, 5, 6, 7, 8, 9, 11, 14))

random_forest <- train(length ~ .,
             data = training,
             method = "rf",
             preProc = c("center", "scale"),
             metric = "RMSE",
             trControl = ctrl,
             tuneGrid = rf_grid,
             ntree = 10,
             importance=TRUE)

rf_metrics = evaluate_model(random_forest, testing)

## Random Forest 
## 
## 36810 samples
##    14 predictor
## 
## Pre-processing: centered (14), scaled (14) 
## Resampling: Cross-Validated (5 fold, repeated 3 times) 
## Summary of sample sizes: 29447, 29449, 29449, 29447, 29448, 29448, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE      
##    2    1.168083  0.5005193  0.9238434
##    3    1.157629  0.5093942  0.9146492
##    4    1.151832  0.5142911  0.9088003
##    5    1.147856  0.5176484  0.9049700
##    6    1.146889  0.5185441  0.9034054
##    7    1.146685  0.5188028  0.9028100
##    8    1.148130  0.5177382  0.9036876
##    9    1.148550  0.5174491  0.9039509
##   11    1.150159  0.5162438  0.9051027
##   14    1.150568  0.5159722  0.9052872
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.
## rf variable importance
## 
##         Overall
## xmonth 100.0000
## chel1   32.0915
## cfin2   23.6747
## chel2   17.8241
## sst     12.8506
## lcop2   12.6688
## totaln  11.1845
## lcop1    9.1186
## cfin1    8.8539
## cumf     7.9361
## fbar     3.9435
## recr     2.4297
## sal      0.2044
## nao      0.0000

## [1] "RMSE na zbiorze testowym: 1.1447769654101 R-Squared na zbiorze testowym 0.520876878600763"

Gradient Boosting Regressor

Podane parametry zostały znaliezone wcześniej wykorzystując grid search, jednakż trwało to bardzo długo i żeby nie ponawiać tych obliczeń parametry zostały tu wpisane na sztywno.

modelLookup(model='gbm')

##   model         parameter                   label forReg forClass
## 1   gbm           n.trees   # Boosting Iterations   TRUE     TRUE
## 2   gbm interaction.depth          Max Tree Depth   TRUE     TRUE
## 3   gbm         shrinkage               Shrinkage   TRUE     TRUE
## 4   gbm    n.minobsinnode Min. Terminal Node Size   TRUE     TRUE
##   probModel
## 1      TRUE
## 2      TRUE
## 3      TRUE
## 4      TRUE

gbm_grid <- expand.grid(
  n.trees=c(10,20,50,100,500,1000),
  shrinkage=c(0.01,0.05,0.1,0.5),
  n.minobsinnode = c(3,5,10),
  interaction.depth=c(1,5,10)
  )

set.seed(23)

params <- data.frame( shrinkage = 0.1, n.trees = 1000, interaction.depth = 10, n.minobsinnode = 10)


gbm <- train(length ~ ., 
             data = training, 
             method = "gbm", 
             metric = "RMSE",
             tuneGrid = params,
             # trControl = ctrl,
             # tuneGrid = gbm_grid,
             preProc = c("center", "scale"),
             verbose = FALSE)

gbm_metrics = evaluate_model(gbm, testing)

## Stochastic Gradient Boosting 
## 
## 36810 samples
##    14 predictor
## 
## Pre-processing: centered (14), scaled (14) 
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 36810, 36810, 36810, 36810, 36810, 36810, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE      
##   1.148978  0.5181491  0.9045936
## 
## Tuning parameter 'n.trees' was held constant at a value of 1000
##  10
## Tuning parameter 'shrinkage' was held constant at a value of
##  0.1
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## gbm variable importance
## 
##         Overall
## sst    100.0000
## recr    27.6999
## xmonth  19.9096
## lcop1   14.6281
## totaln   8.4273
## cfin2    7.0833
## lcop2    6.6431
## chel2    4.8450
## cfin1    3.0231
## sal      0.8983
## nao      0.7995
## cumf     0.6333
## fbar     0.5033
## chel1    0.0000

## [1] "RMSE na zbiorze testowym: 1.13672654284413 R-Squared na zbiorze testowym 0.527527267988041"

ggplot(varImp(gbm))

Podsumowanie regresorów

Jak można zauważyć poniżej najlepszymi modelami okazały się modele oparte na algorytmach Gradient Boosting i Random Forest. Jest nieznaczna różnica pomiędzy tymi algorytmami. Nieco gorzej sprawdza się algorytm random forest. Najgorzej i poniżej oczekiwań działa algorytm Elastic Net, co może być związane z wykorzystaniem regresji liniowej wewnątrz algorytmu.

metrics <- data.frame(
   regressor = c("Elastic Net","Random Forest","Gradient Boosting Regressor"),
   rmse = c(glmnet_metrics[1], rf_metrics[1], gbm_metrics[1]), 
   rsquared = c(glmnet_metrics[2], rf_metrics[2], gbm_metrics[2])
)

metrics

##                     regressor     rmse  rsquared
## 1                 Elastic Net 1.359449 0.3242935
## 2               Random Forest 1.144777 0.5208769
## 3 Gradient Boosting Regressor 1.136727 0.5275273

Analizując istotność danych dla Elastic Net trzema najistotniejszymi atrybutami są fbar, cumf i sst, a najmniej istotnymi chel2, xmonth i sal. Dla algorytmu random forest najistotniejszymi były xmonth, chel1 i cfin2, a najmniej istotne cumf, sal i mao. Najważniejszymi atrybutami dla algorytmu Gradient Boosting były sst, recr i xmonth, a najmniej istotnymi chel1, fbar i sal (wszystkie bardzo niskie).

Można się spodziewać, że atrybut xmonth jest skorelowany z długością (jest to najistotniejszy atrybut dla Random Forest i trzeci najbardziej istotny dla Gradient Boosting). Jednakże nie oznacza to, że jest to wynik wieloletniego trendu, a może na przykład wynikać z wahań związanych z okresami godowymi lub okresami zwiększonego połowu.
Można zauważyć dużą istotność atrybutu sst - najistotniejszy atrybut dla Gradient Boosting, trzeci najbardziej istotny dla Elastic Net. Na wykresie poniżej można zauważyć że średnia długość śledzia jest przeciwnie skorelowana z temperaturą wody przy powierzchni. Może to świadczyć, że jest to powód zmian długości śledzi na przestrzeni lat.

ggplot(aggregated_data, aes(x = time_period)) +
  geom_point(aes(y = sst * 2, colour = "sst")) +
  geom_point(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~.*0.5, name = columns_names["lcop1"])) + 
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) + 
  ggtitle("Długość śledzia wraz z temperaturą pod powierzchnią wody na przestrzeni lat")

Mimo występowania korelacji pomiędzy długością, a oscylacją północnoatlantycką, żaden z regresorów nie uznawał tego atrybutu za istotny. Mniejszą zależność widać w przypadku lcop2 (i chel2).

ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
  geom_line(aes(y = nao + 25, colour = "nao")) +
  geom_line(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~.-25, name = columns_names["nao"])) +
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) +
  ggtitle("Długość śledzia wraz z oscylacja północnoatlantycką na przestrzeni lat")

Często istotnym atrybutem jest atrybut lcop1 (jak i chel1 poprzez korelacje) czyli dostępność planktonu [zagęszczenie widłonogów gat. 1]. Można zauważyć podobny trend zmian tych wartości do trendu zmian wartości długości. Dodatkowo jest on skorelowany z oscylacją północnoatlantycką, przez co ten drugi atrybut może nie występować w istotnych atrybutach.

ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
  geom_line(aes(y = lcop1 / 10 + 25 , colour = "lcop1")) +
  geom_line(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~ (.- 25) * 10, name = columns_names["lcop1"])) +
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) +
  ggtitle("Długość śledzia wraz z zagęszczeniem widłogonów gat. 1 na przestrzeni lat")

ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
  geom_line(aes(y = chel1 / 10 + 25 , colour = "chel1")) +
  geom_line(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~ (.- 25) * 10, name = columns_names["chel1"])) +
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) +
  ggtitle("Długość śledzia wraz z atrubutem chel1 na przestrzeni lat")

ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
  geom_line(aes(y = lcop2 / 10 + 22 , colour = "lcop2")) +
  geom_line(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~ (.- 22) * 10, name = columns_names["lcop2"])) +
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) +
  ggtitle("Długość śledzia wraz z atrubutem lcop2 na przestrzeni lat")

Wartości cumf i rect wydają się nie mieć wpływu na zmianę długości śledzia na przestrzeni lat.

ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
  geom_line(aes(y = cumf * 100 , colour = "cumf")) +
  geom_line(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~. / 100, name = columns_names["cumf"])) +
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) +
  ggtitle("Długość śledzia wraz z rocznym natęzeniem połowów na przestrzeni lat")

ggplot(aggregate_data(completed_data, 1000), aes(x = time_period)) +
  geom_line(aes(y = recr * 0.00005 , colour = "recr")) +
  geom_line(aes(y = length, colour = "length")) +
  scale_y_continuous(sec.axis = sec_axis(~. / 0.00005, name = columns_names["recr"])) +
  scale_colour_manual(values = c("blue", "red")) +
  labs(y = columns_names["length"], x = "Przedział czasowy", colour = "Wartości") +
  theme(legend.position = c(0.25, 0.1)) +
  ggtitle("Długość śledzia wraz z rocznym narybkiem na przestrzeni lat")

Projekt z analizy danych - analiza długości śledzia oceanicznego wyławianego w Europie.

Tomasz Sienkiewicz, Jakub Zdanowski

05 December, 2019